#################################################################################################
# Replication code for
#
# The importance of citizenship for deserving COVID-19 treatment 
# by Rahsaan Maxwell, Marc Helbling, Simon Munzert, and Richard Traunmüller
# 
# Accepted for publication in: Humanities and Social Sciences Communications
# 
# 3.8.2022
# Contact: traunmueller@uni-mannheim.de
#################################################################################################

#################################################################################################
# Load required packages

library(readstata13)
library(dplyr)
library(multiwayvcov)
library(lmtest)
library(RColorBrewer)
library(xtable)
library(arm)

#################################################################################################
# Get data

setwd("~/YOURPATH/")
triage.dat <- read.dta13("TriageItems_w1_to_23.dta",  convert.factors = TRUE, generate.factors = TRUE)

#################################################################################################
# Create Variables

triage.dat$Personenexp <- ifelse(triage.dat$Personenexp==2, 0, 1)

triage.dat$Personenexp[triage.dat$Personenexperiment==1] <- 1  # old version
triage.dat$Personenexp[triage.dat$Personenexperiment==2] <- 0

triage.dat$Personenexp[triage.dat$VAR00001==1] <- 1  # old version
triage.dat$Personenexp[triage.dat$VAR00001==2] <- 0

triage.dat$r.female <- ifelse(triage.dat$v_2==2, 1, 0)
triage.dat$r.age <- (2020-triage.dat$v_3)
triage.dat$r.german <- ifelse(triage.dat$v_9==1, 1, 0)
triage.dat$r.leftright  <- triage.dat$v_121

triage.dat$r.abi  <- ifelse(triage.dat$v_4==4, 1, 0)  
triage.dat$r.child  <- ifelse(triage.dat$v_10==1, 0, 1)
triage.dat$v_311 <- ifelse(triage.dat$v_311==-77, NA, triage.dat$v_311)
triage.dat$r.nation <- triage.dat$v_311
triage.dat$r.immig  <- triage.dat$v_240

triage.dat$Auswahl_Person <- as.numeric(as.factor(triage.dat$Auswahl_Person))-1


# Create Vignettes
triage.1 <- triage.dat[,c("c_0004", "c_0006", "c_0008", "c_0010", "c_0013", "c_0016", "c_0019", "dupl1_v_3", "dupl1_v_3_1_1", "dupl1_v_3_1_2", "dupl1_v_3_1_3", "dupl1_v_3_2_1", "dupl1_v_3_2_2", "dupl1_v_3_2_3", "dupl1_v_3_3_1", "dupl1_v_3_3_2", "dupl1_v_3_3_3", "dupl1_dupl1_v_3_1_1", "dupl1_dupl1_v_3_1_2", "dupl1_dupl1_v_3_1_3", "dupl1_dupl1_v_3_2_1", "dupl1_dupl1_v_3_2_2", "dupl1_dupl1_v_3_2_3", "dupl1_dupl1_v_3_3_1", "dupl1_dupl1_v_3_3_2", "dupl1_dupl1_v_3_3_3", "wave", "Auswahl_Person", "r.female", "r.age", "r.leftright", "r.german", "r.abi", "r.child", "r.nation", "r.immig", "Personenexp")]
triage.2 <- triage.dat[,c("c_0005", "c_0007", "c_0009", "c_0011", "c_0014", "c_0017", "c_0020", "dupl1_v_3", "dupl1_v_3_1_1", "dupl1_v_3_1_2", "dupl1_v_3_1_3", "dupl1_v_3_2_1", "dupl1_v_3_2_2", "dupl1_v_3_2_3", "dupl1_v_3_3_1", "dupl1_v_3_3_2", "dupl1_v_3_3_3", "dupl1_dupl1_v_3_1_1", "dupl1_dupl1_v_3_1_2", "dupl1_dupl1_v_3_1_3", "dupl1_dupl1_v_3_2_1", "dupl1_dupl1_v_3_2_2", "dupl1_dupl1_v_3_2_3", "dupl1_dupl1_v_3_3_1", "dupl1_dupl1_v_3_3_2", "dupl1_dupl1_v_3_3_3", "wave", "Auswahl_Person", "r.female", "r.age", "r.leftright", "r.german", "r.abi", "r.child", "r.nation", "r.immig", "Personenexp")]

triage.1$gender <- as.factor(triage.1$c_0004)
triage.1$age <- as.factor(triage.1$c_0006)
triage.1$job <- as.factor(triage.1$c_0008)
triage.1$chance <- as.factor(triage.1$c_0010)
triage.1$child <- as.factor(triage.1$c_0013)
triage.1$citizen <- as.factor(triage.1$c_0016)
triage.1$crime <- as.factor(triage.1$c_0019)

triage.1$mig2[triage.1$dupl1_v_3_1_1!=-77] <- "German"
triage.1$mig2[triage.1$dupl1_v_3_1_2!=-77] <- "German"
triage.1$mig2[triage.1$dupl1_v_3_1_3!=-77] <- "German"

triage.1$mig2[triage.1$dupl1_v_3_2_1!=-77] <- "EU"
triage.1$mig2[triage.1$dupl1_v_3_2_2!=-77] <- "EU"
triage.1$mig2[triage.1$dupl1_v_3_2_3!=-77] <- "EU"

triage.1$mig2[triage.1$dupl1_v_3_3_1!=-77] <- "Non_EU"
triage.1$mig2[triage.1$dupl1_v_3_3_2!=-77] <- "Non_EU"
triage.1$mig2[triage.1$dupl1_v_3_3_3!=-77] <- "Non_EU"

triage.1$mig2[triage.1$dupl1_dupl1_v_3_1_1!=-77] <- "German"
triage.1$mig2[triage.1$dupl1_dupl1_v_3_1_2!=-77] <- "German"
triage.1$mig2[triage.1$dupl1_dupl1_v_3_1_3!=-77] <- "German"

triage.1$mig2[triage.1$dupl1_dupl1_v_3_2_1!=-77] <- "EU"
triage.1$mig2[triage.1$dupl1_dupl1_v_3_2_2!=-77] <- "EU"
triage.1$mig2[triage.1$dupl1_dupl1_v_3_2_3!=-77] <- "EU"

triage.1$mig2[triage.1$dupl1_dupl1_v_3_3_1!=-77] <- "Non_EU"
triage.1$mig2[triage.1$dupl1_dupl1_v_3_3_2!=-77] <- "Non_EU"
triage.1$mig2[triage.1$dupl1_dupl1_v_3_3_3!=-77] <- "Non_EU"


triage.2$gender <- as.factor(triage.2$c_0005)
triage.2$age <- as.factor(triage.2$c_0007)
triage.2$job <- as.factor(triage.2$c_0009)
triage.2$chance <- as.factor(triage.2$c_0011)
triage.2$child <- as.factor(triage.2$c_0014)
triage.2$citizen <- as.factor(triage.2$c_0017)
triage.2$crime <- as.factor(triage.2$c_0020)

triage.2$mig2[triage.1$dupl1_v_3_1_1!=-77] <- "German"
triage.2$mig2[triage.1$dupl1_v_3_1_2!=-77] <- "EU"
triage.2$mig2[triage.1$dupl1_v_3_1_3!=-77] <- "Non_EU"

triage.2$mig2[triage.1$dupl1_v_3_2_1!=-77] <- "German"
triage.2$mig2[triage.1$dupl1_v_3_2_2!=-77] <- "EU"
triage.2$mig2[triage.1$dupl1_v_3_2_3!=-77] <- "Non_EU"

triage.2$mig2[triage.1$dupl1_v_3_3_1!=-77] <- "German"
triage.2$mig2[triage.1$dupl1_v_3_3_2!=-77] <- "EU"
triage.2$mig2[triage.1$dupl1_v_3_3_3!=-77] <- "Non_EU"

triage.2$mig2[triage.1$dupl1_dupl1_v_3_1_1!=-77] <- "German"
triage.2$mig2[triage.1$dupl1_dupl1_v_3_1_2!=-77] <- "EU"
triage.2$mig2[triage.1$dupl1_dupl1_v_3_1_3!=-77] <- "Non_EU"

triage.2$mig2[triage.1$dupl1_dupl1_v_3_2_1!=-77] <- "German"
triage.2$mig2[triage.1$dupl1_dupl1_v_3_2_2!=-77] <- "EU"
triage.2$mig2[triage.1$dupl1_dupl1_v_3_2_3!=-77] <- "Non_EU"

triage.2$mig2[triage.1$dupl1_dupl1_v_3_3_1!=-77] <- "German"
triage.2$mig2[triage.1$dupl1_dupl1_v_3_3_2!=-77] <- "EU"
triage.2$mig2[triage.1$dupl1_dupl1_v_3_3_3!=-77] <- "Non_EU"


triage.1 <- triage.1[,-c(1:7)]
triage.2 <- triage.2[,-c(1:7)]

triage.1$patient <- 1
triage.2$patient <- 2

triage <- rbind(triage.1, triage.2)

triage$y <- NA
triage$y[triage$dupl1_v_3==1 & triage$patient==1] <- 1
triage$y[triage$dupl1_v_3==2 & triage$patient==2] <- 1
triage$y[triage$Auswahl_Person==1 & triage$patient==1] <- 1
triage$y[triage$Auswahl_Person==2 & triage$patient==2] <- 1
triage$y[triage$dupl1_v_3==1 & triage$patient==2] <- 0
triage$y[triage$dupl1_v_3==2 & triage$patient==1] <- 0
triage$y[triage$Auswahl_Person==1 & triage$patient==2] <- 0
triage$y[triage$Auswahl_Person==2 & triage$patient==1] <- 0

# New Coding Scheme
triage$y[triage$dupl1_v_3_1_1==1 & triage$patient==1] <- 1
triage$y[triage$dupl1_v_3_1_1==2 & triage$patient==2] <- 1
triage$y[triage$dupl1_v_3_1_1==1 & triage$patient==2] <- 0
triage$y[triage$dupl1_v_3_1_1==2 & triage$patient==1] <- 0

triage$y[triage$dupl1_v_3_1_2==1 & triage$patient==1] <- 1
triage$y[triage$dupl1_v_3_1_2==2 & triage$patient==2] <- 1
triage$y[triage$dupl1_v_3_1_2==1 & triage$patient==2] <- 0
triage$y[triage$dupl1_v_3_1_2==2 & triage$patient==1] <- 0

triage$y[triage$dupl1_v_3_1_3==1 & triage$patient==1] <- 1
triage$y[triage$dupl1_v_3_1_3==2 & triage$patient==2] <- 1
triage$y[triage$dupl1_v_3_1_3==1 & triage$patient==2] <- 0
triage$y[triage$dupl1_v_3_1_3==2 & triage$patient==1] <- 0

triage$y[triage$dupl1_v_3_2_1==1 & triage$patient==1] <- 1
triage$y[triage$dupl1_v_3_2_1==2 & triage$patient==2] <- 1
triage$y[triage$dupl1_v_3_2_1==1 & triage$patient==2] <- 0
triage$y[triage$dupl1_v_3_2_1==2 & triage$patient==1] <- 0

triage$y[triage$dupl1_v_3_2_2==1 & triage$patient==1] <- 1
triage$y[triage$dupl1_v_3_2_2==2 & triage$patient==2] <- 1
triage$y[triage$dupl1_v_3_2_2==1 & triage$patient==2] <- 0
triage$y[triage$dupl1_v_3_2_2==2 & triage$patient==1] <- 0

triage$y[triage$dupl1_v_3_2_3==1 & triage$patient==1] <- 1
triage$y[triage$dupl1_v_3_2_3==2 & triage$patient==2] <- 1
triage$y[triage$dupl1_v_3_2_3==1 & triage$patient==2] <- 0
triage$y[triage$dupl1_v_3_2_3==2 & triage$patient==1] <- 0

triage$y[triage$dupl1_v_3_3_1==1 & triage$patient==1] <- 1
triage$y[triage$dupl1_v_3_3_1==2 & triage$patient==2] <- 1
triage$y[triage$dupl1_v_3_3_1==1 & triage$patient==2] <- 0
triage$y[triage$dupl1_v_3_3_1==2 & triage$patient==1] <- 0

triage$y[triage$dupl1_v_3_3_2==1 & triage$patient==1] <- 1
triage$y[triage$dupl1_v_3_3_2==2 & triage$patient==2] <- 1
triage$y[triage$dupl1_v_3_3_2==1 & triage$patient==2] <- 0
triage$y[triage$dupl1_v_3_3_2==2 & triage$patient==1] <- 0

triage$y[triage$dupl1_v_3_3_3==1 & triage$patient==1] <- 1
triage$y[triage$dupl1_v_3_3_3==2 & triage$patient==2] <- 1
triage$y[triage$dupl1_v_3_3_3==1 & triage$patient==2] <- 0
triage$y[triage$dupl1_v_3_2_3==2 & triage$patient==1] <- 0

triage$y[triage$dupl1_dupl1_v_3_1_1==1 & triage$patient==1] <- 1
triage$y[triage$dupl1_dupl1_v_3_1_1==2 & triage$patient==2] <- 1
triage$y[triage$dupl1_dupl1_v_3_1_1==1 & triage$patient==2] <- 0
triage$y[triage$dupl1_dupl1_v_3_1_1==2 & triage$patient==1] <- 0

triage$y[triage$dupl1_dupl1_v_3_1_2==1 & triage$patient==1] <- 1
triage$y[triage$dupl1_dupl1_v_3_1_2==2 & triage$patient==2] <- 1
triage$y[triage$dupl1_dupl1_v_3_1_2==1 & triage$patient==2] <- 0
triage$y[triage$dupl1_dupl1_v_3_1_2==2 & triage$patient==1] <- 0

triage$y[triage$dupl1_dupl1_v_3_1_3==1 & triage$patient==1] <- 1
triage$y[triage$dupl1_dupl1_v_3_1_3==2 & triage$patient==2] <- 1
triage$y[triage$dupl1_dupl1_v_3_1_3==1 & triage$patient==2] <- 0
triage$y[triage$dupl1_dupl1_v_3_1_3==2 & triage$patient==1] <- 0

triage$y[triage$dupl1_dupl1_v_3_2_1==1 & triage$patient==1] <- 1
triage$y[triage$dupl1_dupl1_v_3_2_1==2 & triage$patient==2] <- 1
triage$y[triage$dupl1_dupl1_v_3_2_1==1 & triage$patient==2] <- 0
triage$y[triage$dupl1_dupl1_v_3_2_1==2 & triage$patient==1] <- 0

triage$y[triage$dupl1_dupl1_v_3_2_2==1 & triage$patient==1] <- 1
triage$y[triage$dupl1_dupl1_v_3_2_2==2 & triage$patient==2] <- 1
triage$y[triage$dupl1_dupl1_v_3_2_2==1 & triage$patient==2] <- 0
triage$y[triage$dupl1_dupl1_v_3_2_2==2 & triage$patient==1] <- 0

triage$y[triage$dupl1_dupl1_v_3_2_3==1 & triage$patient==1] <- 1
triage$y[triage$dupl1_dupl1_v_3_2_3==2 & triage$patient==2] <- 1
triage$y[triage$dupl1_dupl1_v_3_2_3==1 & triage$patient==2] <- 0
triage$y[triage$dupl1_dupl1_v_3_2_3==2 & triage$patient==1] <- 0

triage$y[triage$dupl1_dupl1_v_3_3_1==1 & triage$patient==1] <- 1
triage$y[triage$dupl1_dupl1_v_3_3_1==2 & triage$patient==2] <- 1
triage$y[triage$dupl1_dupl1_v_3_3_1==1 & triage$patient==2] <- 0
triage$y[triage$dupl1_dupl1_v_3_3_1==2 & triage$patient==1] <- 0

triage$y[triage$dupl1_dupl1_v_3_3_2==1 & triage$patient==1] <- 1
triage$y[triage$dupl1_dupl1_v_3_3_2==2 & triage$patient==2] <- 1
triage$y[triage$dupl1_dupl1_v_3_3_2==1 & triage$patient==2] <- 0
triage$y[triage$dupl1_dupl1_v_3_3_2==2 & triage$patient==1] <- 0

triage$y[triage$dupl1_dupl1_v_3_3_3==1 & triage$patient==1] <- 1
triage$y[triage$dupl1_dupl1_v_3_3_3==2 & triage$patient==2] <- 1
triage$y[triage$dupl1_dupl1_v_3_3_3==1 & triage$patient==2] <- 0
triage$y[triage$dupl1_dupl1_v_3_2_3==2 & triage$patient==1] <- 0

triage <- triage[, - c(1:19)]

triage$job[triage$job=="Arzt/<c4>rztin"] <- "Arzt/Ärztin"
triage$job[triage$job=="Koch/K<f6>chin"] <- "Koch/Köchin"
triage$job <- droplevels(triage$job)

triage$citizen[triage$wave>=14 & triage$mig2!="German"] <- "ZuwanderIn  mit Aufenthaltsbewilligung"
triage$citizen[triage$wave>=14 & triage$mig2=="German"] <- "Deutscher StaatsbürgerIn"
triage$citizen <- droplevels(triage$citizen)

#############################################################################################
# Descriptive Statistics

des.tab <- rbind(round(c(mean(triage.dat$r.female), sd(triage.dat$r.female), summary(triage.dat$r.female)[c(1,  6)]), 2),
round(c(mean(triage.dat$r.age), sd(triage.dat$r.age), summary(triage.dat$r.age)[c(1,  6)]), 2),
round(c(mean(triage.dat$r.abi), sd(triage.dat$r.abi), summary(triage.dat$r.abi)[c(1,  6)]), 2),
round(c(mean(triage.dat$r.german), sd(triage.dat$r.german), summary(triage.dat$r.german)[c(1,  6)]), 2),
round(c(mean(triage.dat$r.child), sd(triage.dat$r.child), summary(triage.dat$r.child)[c(1,  6)]), 2),
round(c(mean(triage.dat$r.leftright), sd(triage.dat$r.leftright), summary(triage.dat$r.leftright)[c(1,  6)]), 2),
round(c(mean(triage.dat$r.immig), sd(triage.dat$r.immig), summary(triage.dat$r.immig)[c(1,  6)]), 2))

rownames(des.tab) <- c("Female", "Age", "Education (Abitur)", "German", "Has children", "Left-Right-Ideology", "Immigration Attitudes")

des.tab.a <- rbind(round(c(mean(triage.dat$r.female[triage.dat$wave<14], na.rm=T), sd(triage.dat$r.female[triage.dat$wave<14], na.rm=T), summary(triage.dat$r.female[triage.dat$wave<14])[c(1, 6)]), 2),
                 round(c(mean(triage.dat$r.age[triage.dat$wave<14], na.rm=T), sd(triage.dat$r.age[triage.dat$wave<14], na.rm=T), summary(triage.dat$r.age[triage.dat$wave<14])[c(1, 6)]), 2),
                 round(c(mean(triage.dat$r.abi[triage.dat$wave<14], na.rm=T), sd(triage.dat$r.abi[triage.dat$wave<14], na.rm=T), summary(triage.dat$r.abi[triage.dat$wave<14])[c(1, 6)]), 2),
                 round(c(mean(triage.dat$r.german[triage.dat$wave<14], na.rm=T), sd(triage.dat$r.german[triage.dat$wave<14], na.rm=T), summary(triage.dat$r.german[triage.dat$wave<14])[c(1, 6)]), 2),
                 round(c(mean(triage.dat$r.child[triage.dat$wave<14], na.rm=T), sd(triage.dat$r.child[triage.dat$wave<14], na.rm=T), summary(triage.dat$r.child[triage.dat$wave<14])[c(1, 6)]), 2),
                 round(c(mean(triage.dat$r.leftright[triage.dat$wave<14], na.rm=T), sd(triage.dat$r.leftright[triage.dat$wave<14], na.rm=T), summary(triage.dat$r.leftright[triage.dat$wave<14])[c(1, 6)]), 2),
                 round(c(mean(triage.dat$r.immig[triage.dat$wave<14], na.rm=T), sd(triage.dat$r.immig[triage.dat$wave<14], na.rm=T), summary(triage.dat$r.immig[triage.dat$wave<14])[c(1, 6)]), 2))

rownames(des.tab.a) <- c("Female", "Age", "Education (Abitur)", "German", "Has children", "Left-Right-Ideology", "Immigration Attitudes")

des.tab.b <- rbind(round(c(mean(triage.dat$r.female[triage.dat$wave>=14], na.rm=T), sd(triage.dat$r.female[triage.dat$wave>=14], na.rm=T), summary(triage.dat$r.female[triage.dat$wave>=14])[c(1, 6)]), 2),
                   round(c(mean(triage.dat$r.age[triage.dat$wave>=14], na.rm=T), sd(triage.dat$r.age[triage.dat$wave>=14], na.rm=T), summary(triage.dat$r.age[triage.dat$wave>=14])[c(1, 6)]), 2),
                   round(c(mean(triage.dat$r.abi[triage.dat$wave>=14], na.rm=T), sd(triage.dat$r.abi[triage.dat$wave>=14], na.rm=T), summary(triage.dat$r.abi[triage.dat$wave>=14])[c(1, 6)]), 2),
                   round(c(mean(triage.dat$r.german[triage.dat$wave>=14], na.rm=T), sd(triage.dat$r.german[triage.dat$wave>=14], na.rm=T), summary(triage.dat$r.german[triage.dat$wave>=14])[c(1, 6)]), 2),
                   round(c(mean(triage.dat$r.child[triage.dat$wave>=14], na.rm=T), sd(triage.dat$r.child[triage.dat$wave>=14], na.rm=T), summary(triage.dat$r.child[triage.dat$wave>=14])[c(1, 6)]), 2),
                   round(c(mean(triage.dat$r.leftright[triage.dat$wave>=14], na.rm=T), sd(triage.dat$r.leftright[triage.dat$wave>=14], na.rm=T), summary(triage.dat$r.leftright[triage.dat$wave>=14])[c(1, 6)]), 2),
                   round(c(mean(triage.dat$r.immig[triage.dat$wave>=14], na.rm=T), sd(triage.dat$r.immig[triage.dat$wave>=14], na.rm=T), summary(triage.dat$r.immig[triage.dat$wave>=14])[c(1, 6)]), 2))

rownames(des.tab.b) <- c("Female", "Age", "Education (Abitur)", "German", "Has children", "Left-Right-Ideology", "Immigration Attitudes")


xtable(cbind(des.tab, des.tab.a, des.tab.b))

triage.dat$wav <- ifelse(triage.dat$wave>=14, 1, 0)
summary(lm(wav ~ r.female + r.age + r.abi + r.child + r.leftright + r.german + r.immig, data=triage.dat))


####################################################################################################
# Main Effects on Triage Decisions

triage$chance20 <- ifelse(triage$chance=="20 Prozent", 1, 0)
triage$chance50 <- ifelse(triage$chance=="50 Prozent", 1, 0)
triage$job.k <- ifelse(triage$job=="Koch/Köchin", 1, 0)
triage$job.p <- ifelse(triage$job=="PflegerIn", 1, 0)
triage$job.pro <- ifelse(triage$job=="ProfessorIn", 1, 0)
triage$nochild <- ifelse(triage$child=="hat keine schulpflichtige Kinder", 1, 0)

# Run Model
m.0 <- lm(y ~ chance50 + chance20 + age + job.p + job.k + job.pro + crime + gender + nochild + citizen + as.factor(wave), data=triage)
vcov_wave <- cluster.vcov(m.0, cbind(triage$wave))
coeftest(m.0 , vcov_wave)


# Plot Results

labels <- c("Chance of survival 80%", "Chance of survival 50%", "Chance of survival 20%", NA,
            "30 years old", "44 years old", "64 years old",  NA,
            "Medical Doctor", "Nurse", "Cook", "Professor", NA,
            "No criminal record", "Criminal record", NA,
            "Female", "Male", NA,
            "Has school children", "Has no school children", NA,
            "German citizenship", "Immigrant with residence permit")

my.col <- brewer.pal(7, "Spectral")
my.col <- c(rep(my.col[7], 3), NA, rep(my.col[6], 3), NA, rep(my.col[5], 4), NA, rep(my.col[4], 2), NA, rep(my.col[3], 2), NA,
            rep(my.col[2], 2), NA, rep(my.col[1], 2))

coefs <- coeftest(m.0 , vcov_wave)[,1]
coefs <- coefs[-c(1,13:34)]
ci.lo <- coeftest(m.0 , vcov_wave)[,1] - 1.96 * coeftest(m.0 , vcov_wave)[,2]
ci.hi <- coeftest(m.0 , vcov_wave)[,1] + 1.96 * coeftest(m.0 , vcov_wave)[,2]
ci.lo <- ci.lo[-c(1,13:34)]
ci.hi <- ci.hi[-c(1,13:34)]

coefs <- c(0, coefs[1:2], NA, 0, coefs[3:4], NA, 0, coefs[5:7], NA, 0, coefs[8], NA, 0, coefs[9], NA, 0, coefs[10], NA, 0, coefs[11])
ci.lo <- c(0, ci.lo[1:2], NA, 0, ci.lo[3:4], NA, 0, ci.lo[5:7], NA, 0, ci.lo[8], NA, 0, ci.lo[9], NA, 0, ci.lo[10], NA, 0, ci.lo[11])
ci.hi <- c(0, ci.hi[1:2], NA, 0, ci.hi[3:4], NA, 0, ci.hi[5:7], NA, 0, ci.hi[8], NA, 0, ci.hi[9], NA, 0, ci.hi[10], NA, 0, ci.hi[11])

coefs <- rev(coefs)
ci.lo <- rev(ci.lo)
ci.hi <- rev(ci.hi)
labels <- rev(labels)
my.col <- rev(my.col)

par(mar=c(2.5, 8, 1, 2), las=1, mgp=c(1.5, .3, 0), oma=c(2,2,2,1), cex.axis=.8, cex.lab=.8)
plot(coefs, 1:length(coefs), pch=21, xlim=c(-.3, .3), ylim=c(0.2, 24.7), axes=F, cex=.8, ylab="",  xlab="Change in choice probability", bg=my.col, col="black")
grid(lty=1, lwd=.5, col="lightgrey")
axis(2, at=1:length(coefs), labels=labels, col="white")
axis(1, at=seq(-.3, .3, by=.1), labels=c("-30%", "-20%", "-10%", "+/-0%", "+10%", "+20%", "+30%"), col="white")
axis(3, at=seq(-.3, .3, by=.1), labels=c("-30%", "-20%", "-10%", "+/-0%", "+10%", "+20%", "+30%"), col="white")
segments(ci.lo, 1:length(coefs), ci.hi, 1:length(coefs), col="black")
abline(v=0, lty=1, col="grey")
points(coefs, 1:length(coefs), pch=21, bg=my.col, col="black")
arrows(0.02, 24.5, .3, 24.5, length=.05, code=2, col=my.col[24])
arrows(-0.02, 24.5, -.3, 24.5, length=.05, code=2, col=my.col[1])
arrows(0.02, 0.5, .3, 0.5, length=.05, code=2, col=my.col[24])
arrows(-0.02, 0.5, -.3, 0.5, length=.05, code=2, col=my.col[1])
text(0, 24.8, "Higher probability", col=my.col[24], pos=4, cex=.6)
text(0, 24.8, "Lower probability", col=my.col[1], pos=2, cex=.6)


###########################################################################################
# Effect Changes Over Time
# Figure A1

coef.mat <- matrix(NA, 12, 23)
se.mat <- matrix(NA, 12, 23)

for(i in 1:23){
 m <- lm(y ~ chance + gender + age + job + child + citizen + crime, data=triage[triage$wave==i,])
 coef.mat[,i] <- coef(m)
 se.mat[,i] <- se.coef(m)
}

rownames(coef.mat) <- names(coef(m))

par(mfrow=c(3,4), mar=c(2,3,1,1), tck=-.02, mgp=c(1.5, .7, 0))
date <- c(1:7, c(8, 10, 12, 14, 16, 18), c(27, 29, 31, 33, 35, 37, 39, 41, 43, 45))

# "chance50 Prozent"                             
plot(date, coef.mat[2,], ylim=c(-.4, .4), pch=19, main="Chance of survival 50%", xlab="", ylab="Effect on Choice Prob.", xlim=c(-1, 47), axes=F)
grid(lwd=.5, col="lightgrey", lty=1)
segments(date, coef.mat[2,]-1.96*se.mat[2,], date, coef.mat[2,]+1.96*se.mat[2,])
abline(h=0, lty=1, col="lightgrey")
axis(1, at=seq(-1, 47, by=4), labels=c("APR\n2020", "MAY\n2020", "JUN\n2020", "JUL\n2020", "AUG\n2020", "SEP\n2020", "OCT\n2020", "NOV\n2020", "DEC\n2020", "JAN\n2021", "FEB\n2021", "MAR\n2021", "APR\n2021"))
axis(2)
box()


#  "chance80 Prozent"  
plot(date, coef.mat[3,], ylim=c(-.4, .4), pch=19, main="Chance of survival 80%", xlab="", ylab="Effect on Choice Prob.", xlim=c(-1, 47), axes=F)
grid(lwd=.5, col="lightgrey", lty=1)
segments(date, coef.mat[3,]-1.96*se.mat[3,],  date, coef.mat[3,]+1.96*se.mat[3,])
abline(h=0, lty=1, col="lightgrey")
axis(1, at=seq(-1, 47, by=4), labels=c("APR\n2020", "MAY\n2020", "JUN\n2020", "JUL\n2020", "AUG\n2020", "SEP\n2020", "OCT\n2020", "NOV\n2020", "DEC\n2020", "JAN\n2021", "FEB\n2021", "MAR\n2021", "APR\n2021"))
axis(2)
box()


#  "childhat schulpflichtige Kinder"     
plot(date, coef.mat[10,], ylim=c(-.4, .4), pch=19, main="Has school children", xlab="", ylab="Effect on Choice Prob.", xlim=c(-1, 47), axes=F)
grid(lwd=.5, col="lightgrey", lty=1)
segments(date, coef.mat[10,]-1.96*se.mat[10,],  date, coef.mat[10,]+1.96*se.mat[10,])
abline(h=0, lty=1, col="lightgrey")
axis(1, at=seq(-1, 47, by=4), labels=c("APR\n2020", "MAY\n2020", "JUN\n2020", "JUL\n2020", "AUG\n2020", "SEP\n2020", "OCT\n2020", "NOV\n2020", "DEC\n2020", "JAN\n2021", "FEB\n2021", "MAR\n2021", "APR\n2021"))
axis(2)
box()


#  "jobPflegerIn"  
plot(date, coef.mat[8,], ylim=c(-.4, .4), pch=19, main="Nurse", xlab="", ylab="Effect on Choice Prob.", xlim=c(-1, 47), axes=F)
grid(lwd=.5, col="lightgrey", lty=1)
segments(date, coef.mat[8,]-1.96*se.mat[8,],  date, coef.mat[8,]+1.96*se.mat[8,])
abline(h=0, lty=1, col="lightgrey")
axis(1, at=seq(-1, 47, by=4), labels=c("APR\n2020", "MAY\n2020", "JUN\n2020", "JUL\n2020", "AUG\n2020", "SEP\n2020", "OCT\n2020", "NOV\n2020", "DEC\n2020", "JAN\n2021", "FEB\n2021", "MAR\n2021", "APR\n2021"))
axis(2)
box()


#  "jobKoch/Köchin"       
plot(date, coef.mat[7,], ylim=c(-.4, .4), pch=19, main="Cook", xlab="", ylab="Effect on Choice Prob.", xlim=c(-1, 47), axes=F)
grid(lwd=.5, col="lightgrey", lty=1)
segments(date, coef.mat[7,]-1.96*se.mat[7,],  date, coef.mat[7,]+1.96*se.mat[7,])
abline(h=0, lty=1, col="lightgrey")
axis(1, at=seq(-1, 47, by=4), labels=c("APR\n2020", "MAY\n2020", "JUN\n2020", "JUL\n2020", "AUG\n2020", "SEP\n2020", "OCT\n2020", "NOV\n2020", "DEC\n2020", "JAN\n2021", "FEB\n2021", "MAR\n2021", "APR\n2021"))
axis(2)
box()


#  "jobProfessorIn"   
plot(date, coef.mat[9,], ylim=c(-.4, .4), pch=19, main="Professor", xlab="", ylab="Effect on Choice Prob.", xlim=c(-1, 47), axes=F)
grid(lwd=.5, col="lightgrey", lty=1)
segments(date, coef.mat[9,]-1.96*se.mat[9,],  date, coef.mat[9,]+1.96*se.mat[9,])
abline(h=0, lty=1, col="lightgrey")
axis(1, at=seq(-1, 47, by=4), labels=c("APR\n2020", "MAY\n2020", "JUN\n2020", "JUL\n2020", "AUG\n2020", "SEP\n2020", "OCT\n2020", "NOV\n2020", "DEC\n2020", "JAN\n2021", "FEB\n2021", "MAR\n2021", "APR\n2021"))
axis(2)
box()


#  "age44"  
plot(date, coef.mat[5,], ylim=c(-.4, .4), pch=19, main="44 years old", xlab="", ylab="Effect on Choice Prob.", xlim=c(-1, 47), axes=F)
grid(lwd=.5, col="lightgrey", lty=1)
segments(date, coef.mat[5,]-1.96*se.mat[5,],  date, coef.mat[5,]+1.96*se.mat[5,])
abline(h=0, lty=1, col="lightgrey")
axis(1, at=seq(-1, 47, by=4), labels=c("APR\n2020", "MAY\n2020", "JUN\n2020", "JUL\n2020", "AUG\n2020", "SEP\n2020", "OCT\n2020", "NOV\n2020", "DEC\n2020", "JAN\n2021", "FEB\n2021", "MAR\n2021", "APR\n2021"))
axis(2)
box()


#  "age64" 
plot(date, coef.mat[6,], ylim=c(-.4, .4), pch=19, main="64 years old", xlab="", ylab="Effect on Choice Prob.", xlim=c(-1, 47), axes=F)
grid(lwd=.5, col="lightgrey", lty=1)
segments(date, coef.mat[6,]-1.96*se.mat[6,],  date, coef.mat[6,]+1.96*se.mat[6,])
abline(h=0, lty=1, col="lightgrey")
axis(1, at=seq(-1, 47, by=4), labels=c("APR\n2020", "MAY\n2020", "JUN\n2020", "JUL\n2020", "AUG\n2020", "SEP\n2020", "OCT\n2020", "NOV\n2020", "DEC\n2020", "JAN\n2021", "FEB\n2021", "MAR\n2021", "APR\n2021"))
axis(2)
box()


#  "genderMann" 
plot(date, coef.mat[4,], ylim=c(-.4, .4), pch=19, main="Male", xlab="", ylab="Effect on Choice Prob.", xlim=c(-1, 47), axes=F)
grid(lwd=.5, col="lightgrey", lty=1)
segments(date, coef.mat[4,]-1.96*se.mat[4,],  date, coef.mat[4,]+1.96*se.mat[4,])
abline(h=0, lty=1, col="lightgrey")
axis(1, at=seq(-1, 47, by=4), labels=c("APR\n2020", "MAY\n2020", "JUN\n2020", "JUL\n2020", "AUG\n2020", "SEP\n2020", "OCT\n2020", "NOV\n2020", "DEC\n2020", "JAN\n2021", "FEB\n2021", "MAR\n2021", "APR\n2021"))
axis(2)
box()


#  "citizenZuwanderIn  mit Aufenthaltsbewilligung"
plot(date, coef.mat[11,], ylim=c(-.4, .4), pch=19, main="Immigrant", xlab="", ylab="Effect on Choice Prob.", xlim=c(-1, 47), axes=F)
grid(lwd=.5, col="lightgrey", lty=1)
segments(date, coef.mat[11,]-1.96*se.mat[11,],  date, coef.mat[11,]+1.96*se.mat[11,])
abline(h=0, lty=1, col="lightgrey")
axis(1, at=seq(-1, 47, by=4), labels=c("APR\n2020", "MAY\n2020", "JUN\n2020", "JUL\n2020", "AUG\n2020", "SEP\n2020", "OCT\n2020", "NOV\n2020", "DEC\n2020", "JAN\n2021", "FEB\n2021", "MAR\n2021", "APR\n2021"))
axis(2)
box()

# Test for time effect
summary(lm(coef.mat[11,] ~ c(1:23) + time))
time <- c(rep(0, 13), rep(1, 10))
summary(lm(coef.mat[11,] ~ time))

#  "crimevorbestraft"   
plot(date, coef.mat[12,], ylim=c(-.4, .4), pch=19, main="Criminal record", xlab="", ylab="Effect on Choice Prob.", xlim=c(-1, 47), axes=F)
grid(lwd=.5, col="lightgrey", lty=1)
segments(date, coef.mat[12,]-1.96*se.mat[12,],  date, coef.mat[12,]+1.96*se.mat[12,])
abline(h=0, lty=1, col="lightgrey")
axis(1, at=seq(-1, 47, by=4), labels=c("APR\n2020", "MAY\n2020", "JUN\n2020", "JUL\n2020", "AUG\n2020", "SEP\n2020", "OCT\n2020", "NOV\n2020", "DEC\n2020", "JAN\n2021", "FEB\n2021", "MAR\n2021", "APR\n2021"))
axis(2)
box()

###########################################################################################
# Variation Across Respondents

m.9a <- lm(y ~ chance + gender + age + job + child + citizen + crime + age*r.age, data=triage)
vcov_wave <- cluster.vcov(m.9a, triage$wave)
coeftest(m.9a , vcov_wave)

m.9b <- lm(y ~ chance + gender + age + job + child + citizen + crime + gender*r.female, data=triage)
vcov_wave <- cluster.vcov(m.9b, triage$wave)
coeftest(m.9b , vcov_wave)

m.9c <- lm(y ~ chance + gender + age + job + child + citizen + crime + job*r.abi, data=triage)
vcov_wave <- cluster.vcov(m.9c, triage$wave)
coeftest(m.9c, vcov_wave)

m.9c <- lm(y ~ chance + gender + age + job + child + citizen + crime + job*r.abi, data=triage)
vcov_wave <- cluster.vcov(m.9c, triage$wave)
coeftest(m.9c, vcov_wave)

m.9d <- lm(y ~ chance + gender + age + job + child + citizen + crime + child*r.child, data=triage)
vcov_wave <- cluster.vcov(m.9d, triage$wave)
coeftest(m.9d, vcov_wave)

m.9e <- lm(y ~ chance + gender + age + job + child + citizen + crime + citizen*r.german, data=triage)
vcov_wave <- cluster.vcov(m.9e, triage$wave)
coeftest(m.9e, vcov_wave)


#############################################################################################
# Figure A4

sim.9a <- sim(m.9a, 1000)

my.col <- brewer.pal(3, "Set2")
par(mfrow=c(2,2))
par(mar=c(3,3,2,1),  mgp=c(2, .5, 0), cex.axis=.8, cex.lab=.8)
curve(coef(m.9a)["age44"] + coef(m.9a)["age44:r.age"]*x, xlim=c(18, 80), ylim=c(-.2, .2), axes=F, xlab="R Age in Years", ylab="Marginal Effect of Age", col=my.col[2], lwd=1.5)
for(i in 1:1000){
  curve(coef(sim.9a)[i,"age44"] + coef(sim.9a)[i, "age44:r.age"]*x, col=my.col[2], lwd=0.5, add=T)
}

grid(lty=1, lwd=.5, col="lightgrey")
axis(1, col="white")
axis(2, col="white")

curve(coef(m.9a)["age64"] + coef(m.9a)["age64:r.age"]*x, add=T, col=my.col[3], lwd=1.5)
for(i in 1:1000){
  curve(coef(sim.9a)[i,"age64"] + coef(sim.9a)[i, "age64:r.age"]*x, col=my.col[3], lwd=0.5, add=T)
}

abline(h=0, lty=2)

#text(8, .09, "Patient: 30 years", cex=.8, pos=3, col=my.col[1])
text(30, .05, "Patient: 44 years", cex=.8, pos=4, col=my.col[2])
text(30, -.15, "Patient: 64 years", cex=.8, pos=4, col=my.col[3])
mtext("a)", side=3, line=1, adj=0, font=2, cex=.8)


coefs.9b <- c(coef(m.9b)["genderMann"], 
              coef(m.9b)["genderMann"] + coef(m.9b)["genderMann:r.female"])

sim.9b <- arm::sim(m.9b, 1000)

plot(c(1:2), coefs.9b, pch=19, col=c(my.col[3], my.col[2]), xlim=c(0.5, 2.5), ylim=c(-.1, .1), ylab="Marginal Effect of Male", xlab="", axes=F)
axis(2, col="white")
axis(1, at=c(1, 2), labels=c("R Female", "R Male"), col="white")
grid(lty=1, lwd=.5, col="lightgrey")
segments(1, quantile(coef(sim.9b)[, "genderMann"], .025), 1, quantile(coef(sim.9b)[, "genderMann"], .975), my.col[3])
segments(2, quantile(coef(sim.9b)[,"genderMann"] + coef(sim.9b)[,"genderMann:r.female"], .025), 2, quantile(coef(sim.9b)[,"genderMann"] +coef(sim.9b)[,"genderMann:r.female"], .975), my.col[2])
#box()
abline(h=0, lty=2)
mtext("b)", side=3, line=1, adj=0, font=2, cex=.8)


coefs.9c <- c(coef(m.9c)["jobKoch/Köchin"], coef(m.9c)["jobPflegerIn"], coef(m.9c)["jobProfessorIn"],
              coef(m.9c)["jobKoch/Köchin"]+coef(m.9c)["jobKoch/Köchin:r.abi"], coef(m.9c)["jobPflegerIn"]+coef(m.9c)["jobPflegerIn:r.abi"], coef(m.9c)["jobProfessorIn"]+coef(m.9c)["jobProfessorIn:r.abi"])
              
sim.9c <- arm::sim(m.9c, 1000)

plot(c(.75, 1, 1.25), coefs.9c[1:3], pch=19, col=c(my.col[3]), xlim=c(0.5, 2.5), ylim=c(-.2, .2), ylab="Marginal Effect of Job (Ref: Doctor)", xlab="", axes=F)
axis(2, col="white")
axis(1, at=c(1:2), labels=c("R No College", "R College"), col="white")
grid(lty=1, lwd=.5, col="lightgrey")
segments(.75, quantile(coef(sim.9c)[, "jobKoch/Köchin"], .025), .75, quantile(coef(sim.9c)[, "jobKoch/Köchin"], .975), my.col[3])
segments(1, quantile(coef(sim.9c)[, "jobPflegerIn"], .025), 1, quantile(coef(sim.9c)[, "jobPflegerIn"], .975), my.col[3])
segments(1.25, quantile(coef(sim.9c)[, "jobProfessorIn"], .025), 1.25, quantile(coef(sim.9c)[, "jobProfessorIn"], .975), my.col[3])

points(c(1.75, 2, 2.25), coefs.9c[4:6], pch=19, col=c(my.col[2]))
segments(1.75, quantile(coef(sim.9c)[, "jobKoch/Köchin"]+coef(sim.9c)[, "jobKoch/Köchin:r.abi"], .025), 1.75, quantile(coef(sim.9c)[, "jobKoch/Köchin"]+coef(sim.9c)[, "jobKoch/Köchin:r.abi"], .975), my.col[2])
segments(2, quantile(coef(sim.9c)[, "jobPflegerIn"]+coef(sim.9c)[, "jobPflegerIn:r.abi"], .025), 2, quantile(coef(sim.9c)[, "jobPflegerIn"]+coef(sim.9c)[, "jobPflegerIn:r.abi"], .975), my.col[2])
segments(2.25, quantile(coef(sim.9c)[, "jobProfessorIn"]+coef(sim.9c)[, "jobProfessorIn:r.abi"], .025), 2.25, quantile(coef(sim.9c)[, "jobProfessorIn"]+coef(sim.9c)[, "jobProfessorIn:r.abi"], .975), my.col[2])

text(c(.75, 1, 1.25), c(.05, .1, .0), c("Cook", "Nurse", "Professor"), col=c(my.col[3]), cex=.6, pos=1)
text(c(1.75, 2, 2.25), c(.01, .05, .02), c("Cook", "Nurse", "Professor"), col=c(my.col[2]), cex=.6, pos=1)
mtext("c)", side=3, line=1, adj=0, font=2, cex=.8)

#box()
abline(h=0, lty=2)


coefs.9d <- c(coef(m.9d)["childhat schulpflichtige Kinder"], 
              coef(m.9d)["childhat schulpflichtige Kinder"] + coef(m.9d)["childhat schulpflichtige Kinder:r.child"])

sim.9d <- arm::sim(m.9d, 1000)

plot(c(1:2), coefs.9d, pch=19, col=c(my.col[3], my.col[2]), xlim=c(0.5, 2.5), ylim=c(.1, .3), ylab="Marginal Effect of Children", xlab="", axes=F)
axis(2, col="white")
axis(1, at=c(1, 2), labels=c("R Has No Children", "R Has Children"), col="white")
grid(lty=1, lwd=.5, col="lightgrey")
segments(1, quantile(coef(sim.9d)[, "childhat schulpflichtige Kinder"], .025), 1, quantile(coef(sim.9d)[, "childhat schulpflichtige Kinder"], .975), my.col[3])
segments(2, quantile(coef(sim.9d)[,"childhat schulpflichtige Kinder"] + coef(sim.9d)[,"childhat schulpflichtige Kinder:r.child"], .025), 2, quantile(coef(sim.9d)[,"childhat schulpflichtige Kinder"] +coef(sim.9d)[,"childhat schulpflichtige Kinder:r.child"], .975), my.col[2])
#box()
abline(h=0, lty=2)
mtext("d)", side=3, line=1, adj=0, font=2, cex=.8)

##############################################################################################
# Figure A2

coefs.9e <- c(coef(m.9e)["citizenZuwanderIn  mit Aufenthaltsbewilligung"], 
              coef(m.9e)["citizenZuwanderIn  mit Aufenthaltsbewilligung"] + coef(m.9e)["citizenZuwanderIn  mit Aufenthaltsbewilligung:r.german"])

sim.9e <- arm::sim(m.9e, 1000)

plot(c(1:2), coefs.9e, pch=19, col=c(my.col[3], my.col[2]), xlim=c(0.5, 2.5), ylim=c(-.2, .2), ylab="Marginal Effect of Immgrant Status", xlab="", axes=F)
axis(2, col="white")
axis(1, at=c(1, 2), labels=c("R Immigrant Resident", "R German Citizen"), col="white")
grid(lty=1, lwd=.5, col="lightgrey")
segments(1, quantile(coef(sim.9e)[, "citizenZuwanderIn  mit Aufenthaltsbewilligung"], .025), 1, quantile(coef(sim.9e)[, "citizenZuwanderIn  mit Aufenthaltsbewilligung"], .975), my.col[3])
segments(2, quantile(coef(sim.9e)[,"citizenZuwanderIn  mit Aufenthaltsbewilligung"] + coef(sim.9e)[,"citizenZuwanderIn  mit Aufenthaltsbewilligung:r.german"], .025), 2, quantile(coef(sim.9e)[,"citizenZuwanderIn  mit Aufenthaltsbewilligung"] + coef(sim.9e)[,"citizenZuwanderIn  mit Aufenthaltsbewilligung:r.german"], .975), my.col[2])
#box()
abline(h=0, lty=2)

##############################################################################################
# Figure 2

m.11a <- lm(y ~ chance + gender + age + job + child + citizen + crime + citizen*r.nation, data=triage)
vcov_wave <- cluster.vcov(m.11a, triage$wave)
coeftest(m.11a, vcov_wave)

m.11b <- lm(y ~ chance + gender + age + job + child + citizen + crime + citizen*r.immig, data=triage)
vcov_wave <- cluster.vcov(m.11b, triage$wave)
coeftest(m.11b, vcov_wave)

sim.11a <- sim(m.11a, 1000)
par(mar=c(3,3,2,1), mfrow=c(2,1))
curve(coef(m.11a)["citizenZuwanderIn  mit Aufenthaltsbewilligung"] + coef(m.11a)["citizenZuwanderIn  mit Aufenthaltsbewilligung:r.nation"]*x, xlim=c(0, 100), ylim=c(-.30, .3),  axes=F, xlab="R Feels Close Toward Germans", ylab="Marginal Effect of Immigrant", col=my.col[2], lwd=1.5)
grid(lty=1, lwd=.5, col="lightgrey")
axis(1, col="white")
axis(2, col="white")
abline(h=0, lty=1, col="grey")
for(i in 1:1000){
  curve(coef(sim.11a)[i,"citizenZuwanderIn  mit Aufenthaltsbewilligung"] + coef(sim.11a)[i,"citizenZuwanderIn  mit Aufenthaltsbewilligung:r.nation"]*x, add=T, col=my.col[2], lwd=1)
}

sim.11b <- sim(m.11b, 1000)
curve(coef(m.11b)["citizenZuwanderIn  mit Aufenthaltsbewilligung"] + coef(m.11b)["citizenZuwanderIn  mit Aufenthaltsbewilligung:r.immig"]*x, xlim=c(0, 8), ylim=c(-.30, .3),  axes=F, xlab="R Critical of Immigration", ylab="Marginal Effect of Immigrant", col=my.col[2], lwd=1.5)
grid(lty=1, lwd=.5, col="lightgrey")
axis(1, col="white")
axis(2, col="white")
abline(h=0, lty=1, col="grey")
for(i in 1:1000){
  curve(coef(sim.11b)[i,"citizenZuwanderIn  mit Aufenthaltsbewilligung"] + coef(sim.11b)[i,"citizenZuwanderIn  mit Aufenthaltsbewilligung:r.immig"]*x, add=T, col=my.col[2], lwd=1)
}

##############################################################################################


m.12 <- lm(y ~ chance + gender + age + job + child + citizen + crime + r.leftright + job:r.leftright + age:r.leftright + crime:r.leftright + child:r.leftright, data=triage)
vcov_wave <- cluster.vcov(m.12, triage$wave)
coeftest(m.12, vcov_wave)

m.12a <- lm(y ~ chance + gender + age + job + child + citizen + crime + job*r.leftright, data=triage)
vcov_wave <- cluster.vcov(m.12a, triage$wave)
coeftest(m.12a, vcov_wave)

m.12b <- lm(y ~ chance + gender + age*r.leftright + job + child + citizen + crime + job, data=triage)
vcov_wave <- cluster.vcov(m.12b, triage$wave)
coeftest(m.12b, vcov_wave)

m.12d <- lm(y ~ chance + gender + age + job + child*r.leftright + citizen + crime + job, data=triage)
vcov_wave <- cluster.vcov(m.12d, triage$wave)
coeftest(m.12d, vcov_wave)

m.12c <- lm(y ~ chance + gender + age + job + child + citizen + crime*r.leftright + job, data=triage)
vcov_wave <- cluster.vcov(m.12c, triage$wave)
coeftest(m.12c, vcov_wave)

m.12e <- lm(y ~ chance + gender + age + job + child + citizen*r.leftright + crime, data=triage)
vcov_wave <- cluster.vcov(m.12e, triage$wave)
coeftest(m.12e, vcov_wave)




coef(m.12a)

m.12b <- lm(y ~ chance + gender + age + job + child + citizen + crime + r.leftright + I(r.leftright^2) + age:r.leftright + age:I(r.leftright^2), data=triage)
vcov_wave <- cluster.vcov(m.12b, triage$wave)
coeftest(m.12b, vcov_wave)


my.col <- brewer.pal(3, "Set2")

sim.12a <- sim(m.12a, 1000)
par(mar=c(3,3,2,1), mfrow=c(2,2))
curve(coef(m.12a)["jobKoch/Köchin"] + coef(m.12a)["jobKoch/Köchin:r.leftright"]*x, xlim=c(0, 10), ylim=c(-.30, .20), axes=F, xlab="Left-Right Ideology", ylab="Marginal Effect of Job (Ref: Doctor)", col=my.col[2], lwd=1.5)
grid(lty=1, lwd=.5, col="lightgrey")
axis(1, col="white")
axis(2, col="white")
abline(h=0, lty=1, col="grey")
curve(coef(m.12a)["jobPflegerIn"] + coef(m.12a)["jobPflegerIn:r.leftright"]*x, add=T, col=my.col[3], lwd=1.5)
curve(coef(m.12a)["jobProfessorIn"] + coef(m.12a)["jobProfessorIn:r.leftright"]*x, add=T, col=my.col[1], lwd=1.5)

for(i in 1:1000){
  curve(coef(sim.12a)[i,"jobPflegerIn"] + coef(sim.12a)[i,"jobPflegerIn:r.leftright"]*x, add=T, col=my.col[3], lwd=1)
}

for(i in 1:1000){
  curve(coef(sim.12a)[i,"jobProfessorIn"] + coef(sim.12a)[i,"jobProfessorIn:r.leftright"]*x, add=T, col=my.col[1], lwd=1)
}

for(i in 1:1000){
  curve(coef(sim.12a)[i,"jobKoch/Köchin"] + coef(sim.12a)[i,"jobKoch/Köchin:r.leftright"]*x, add=T, col=my.col[2], lwd=1)
}

text(8, -.15, "Patient: Cook", cex=.8, pos=, col=my.col[2])
text(8, .06, "Patient: Nurse", cex=.8, pos=3, col=my.col[3])
text(2, -.2, "Patient: Professor", cex=.8, pos=3, col=my.col[1])
mtext("a)", side=3, line=1, adj=0, font=2)


sim.12b <- sim(m.12b, 1000)
par(mar=c(3,3,2,1))
curve(coef(m.12b)["age44"] + coef(m.12b)["age44:r.leftright"]*x, xlim=c(0, 10), ylim=c(-.30, .20), axes=F, xlab="Left-Right Ideology", ylab="Marginal Effect of Age", col=my.col[2], lwd=1.5)
grid(lty=1, lwd=.5, col="lightgrey")
axis(1, col="white")
axis(2, col="white")
abline(h=0, lty=1, col="grey")
curve(coef(m.12b)["age64"] + coef(m.12b)["age64:r.leftright"]*x, add=T, col=my.col[3], lwd=1.5)
#curve(coef(m.12b)["r.leftright"]*x, add=T, col=my.col[1])

for(i in 1:1000){
  curve(coef(sim.12b)[i,"age44"] + coef(sim.12b)[i,"age44:r.leftright"]*x, add=T, col=my.col[2], lwd=1)
}

for(i in 1:1000){
  curve(coef(sim.12b)[i,"age64"] + coef(sim.12b)[i,"age64:r.leftright"]*x, add=T, col=my.col[3], lwd=1)
}


#text(8, .09, "Patient: 30 years", cex=.8, pos=3, col=my.col[1])
text(8, .1, "Patient: 44 years", cex=.8, pos=, col=my.col[2])
text(8, -.21, "Patient: 64 years", cex=.8, pos=3, col=my.col[3])
mtext("b)", side=3, line=1, adj=0, font=2)


#m.12c <- lm(y ~ chance + gender + age + job + child + citizen + crime + r.leftright + crime:r.leftright, data=triage)
#vcov_wave <- cluster.vcov(m.12c, triage$wave)
#coeftest(m.12c, vcov_wave)

sim.12c<- sim(m.12c, 1000)
par(mar=c(3,3,2,1))
curve(coef(m.12c)["crimevorbestraft"] + coef(m.12c)["crimevorbestraft:r.leftright"]*x, xlim=c(0, 10),  ylim=c(-.30, 0), axes=F, xlab="Left-Right Ideology", ylab="Marginal Effect of Criminal Record", col=my.col[2], lwd=1.5)
grid(lty=1, lwd=.5, col="lightgrey")
axis(1, col="white")
axis(2, col="white")
abline(h=0, lty=1, col="grey")

for(i in 1:1000){
  curve(coef(sim.12c)[i,"crimevorbestraft"] + coef(sim.12c)[i,"crimevorbestraft:r.leftright"]*x, add=T, col=my.col[2], lwd=1)
}
mtext("c)", side=3, line=1, adj=0, font=2)

#text(5, -.15, "Patient: convicted", cex=.8, pos=3, col=my.col[2])


#m.12d <- lm(y ~ chance + gender + age + job + child + citizen + crime + r.leftright + child:r.leftright, data=triage)
#vcov_wave <- cluster.vcov(m.12d, triage$wave)
#coeftest(m.12d, vcov_wave)
sim.12d <- sim(m.12d, 1000)
par(mar=c(3,3,2,1))
curve(coef(m.12d)["childhat schulpflichtige Kinder"] + coef(m.12d)["childhat schulpflichtige Kinder:r.leftright"]*x, xlim=c(0, 10),  ylim=c(0, .50), axes=F, xlab="R Left-Right Ideology", ylab="Marginal Effect of Having Children", col=my.col[2], lwd=1.5)
grid(lty=1, lwd=.5, col="lightgrey")
axis(1, col="white")
axis(2, col="white")
abline(h=0, lty=1, col="grey")
for(i in 1:1000){
  curve(coef(sim.12d)[i,"childhat schulpflichtige Kinder"] + coef(sim.12d)[i,"childhat schulpflichtige Kinder:r.leftright"]*x, add=T, col=my.col[2], lwd=1)
}
mtext("d)", side=3, line=1, adj=0, font=2)


text(7, .05, "Patient: has school children", cex=.8, pos=3, col=my.col[2])

m.12e <- lm(y ~ chance + gender + age + job + child + citizen + crime + r.leftright + citizen:r.leftright, data=triage)
vcov_wave <- cluster.vcov(m.12e, triage$wave)
coeftest(m.12e, vcov_wave)

par(mar=c(3,3,2,1))
curve(coef(m.12e)["citizenZuwanderIn  mit Aufenthaltsbewilligung"] + coef(m.12e)["citizenZuwanderIn  mit Aufenthaltsbewilligung:r.leftright"]*x, xlim=c(0, 10),  ylim=c(-.50, .2), axes=F, xlab="Left-Right Ideology", ylab="Marginal Effect: Immigrant", col=my.col[2], lwd=1.5)
grid(lty=1, lwd=.5, col="lightgrey")
axis(1, col="white")
axis(2, col="white")
abline(h=0, lty=1, col="grey")

text(7, -.1, "Patient: migrant", cex=.8, pos=3, col=my.col[2])

m.12f <- lm(y ~ chance + gender + age + job + child + citizen + crime + r.immig + citizen:r.immig, data=triage)
vcov_wave <- cluster.vcov(m.12f, triage$wave)
coeftest(m.12f, vcov_wave)

par(mar=c(3,3,2,1))
curve(coef(m.12f)["citizenZuwanderIn  mit Aufenthaltsbewilligung"] + coef(m.12f)["citizenZuwanderIn  mit Aufenthaltsbewilligung:r.immig"]*x, xlim=c(0, 10),  ylim=c(-.50, .2), axes=F, xlab="Anti Immigration", ylab="Marginal Effect: Immigrant", col=my.col[2], lwd=1.5)
grid(lty=1, lwd=.5, col="lightgrey")
axis(1, col="white")
axis(2, col="white")
abline(h=0, lty=1, col="grey")

##############################################################################################
# Variation Among Patients
# Figure 3

triage$superint <- interaction(triage$chance, triage$gender, triage$age, triage$job, triage$child, triage$crime)
m.superint <- lmer(y ~ citizen + (1 + citizen | superint) + as.factor(wave) , data=triage)
summary(m.superint)

par(yaxs="i", mfrow=c(1,1), mar=c(3,3,2,2), cex.lab=.8, cex.axis=.8)
hist(coef(m.superint)$superint[,2], breaks=20, xlim=c(-.20, -.05), main="", xlab="Immigrant Effects (AMCEs)")


##############################################################################################
# EU vs. Non-EU Immigrants
# Figure 4

m.14.b <- lm(y ~ mig2 + as.factor(wave), data=triage[triage$wave>=14, ])
vcov_wave <- cluster.vcov(m.14.b, cbind(triage$wave[triage$wave>=14]))
coeftest(m.14.b , vcov_wave)

coefs.14b <- c(coef(m.14.b)["(Intercept)"], 
               coef(m.14.b)["(Intercept)"] + coef(m.14.b)["mig2Non_EU"],
               coef(m.14.b)["(Intercept)"] + coef(m.14.b)["mig2German"])
sim.14b <- sim(m.14.b, 1000)

plot(c(1:3), coefs.14b, pch=19, col=c(my.col[3]), xlim=c(0.5, 3.5), ylim=c(.4, .6), ylab="Choice Probability", xlab="", axes=F)
axis(2, col="white")
axis(1, at=c(1:3), labels=c("EU Immigrant", "Non-EU Immigrant", "German Citizen"), col="white")
grid(lty=1, lwd=.5, col="lightgrey")
segments(1, quantile(coef(sim.14b)[, "(Intercept)"], .025), 1, quantile(coef(sim.14b)[, "(Intercept)"], .975), my.col[3])
segments(3, quantile(coef(sim.14b)[,"(Intercept)"] + coef(sim.14b)[,"mig2German"], .025), 3, quantile(coef(sim.14b)[,"(Intercept)"] +coef(sim.14b)[,"mig2German"], .975), my.col[3])
segments(2, quantile(coef(sim.14b)[,"(Intercept)"] + coef(sim.14b)[,"mig2Non_EU"], .025), 2, quantile(coef(sim.14b)[,"(Intercept)"] +coef(sim.14b)[,"mig2Non_EU"], .975), my.col[3])

#############################################################################################
# Test of Different Coding Versions
# Figure A3

# Old Version
m.0 <- lm(y ~ chance + gender + age + job + child + citizen + crime + as.factor(wave), data=triage[triage$wave>14 & triage$Personenexp==1,])
summary(m.0)
vcov_wave <- cluster.vcov(m.0, cbind(triage$wave[triage$wave>14]))
coeftest(m.0 , vcov_wave)

coefs <- coeftest(m.0 , vcov_wave)[,1]
coefs <- coefs[-c(1, 13:20)]
ci.lo <- coeftest(m.0 , vcov_wave)[,1] - 1.96 * coeftest(m.0 , vcov_wave)[,2]
ci.hi <- coeftest(m.0 , vcov_wave)[,1] + 1.96 * coeftest(m.0 , vcov_wave)[,2]
ci.lo <- ci.lo[-c(1, 13:20)]
ci.hi <- ci.hi[-c(1, 13:20)]


labels <- c("Chance of survival 50%", "Chance of survival 80%",
            "Male", "44 years old", "64 years old", "Cook", "Nurse", "Professor", "Has school children", "Immigrant", "Criminal record")

par(mar=c(2.5, 8, 1, 2), las=1, mgp=c(1.5, .3, 0), oma=c(2,2,2,1), cex.axis=.8, cex.lab=.8)
plot(coefs, 1:length(coefs), pch=19, xlim=c(-.3, .3), ylim=c(0.5, 11.5), axes=F, cex=.8, ylab="",  xlab="Change in choice probability", col="red")
grid(lty=1, lwd=.5, col="lightgrey")
axis(2, at=1:length(coefs), labels=labels, col="white")
axis(1, at=seq(-.3, .3, by=.1), labels=c("-30%", "-20%", "-10%", "+/-0%", "+10%", "+20%", "+30%"), col="white")
axis(3, at=seq(-.3, .3, by=.1), labels=c("-30%", "-20%", "-10%", "+/-0%", "+10%", "+20%", "+30%"), col="white")
segments(ci.lo, 1:length(coefs), ci.hi, 1:length(coefs), col="red")
abline(v=0, lty=2)

# New
m.0 <- lm(y ~ chance + gender + age + job + child + citizen + crime + as.factor(wave), data=triage[triage$wave>14 & triage$Personenexp==0,])
summary(m.0)
vcov_wave <- cluster.vcov(m.0, cbind(triage$wave[triage$wave>14 & triage$Personenexp==0]))
coeftest(m.0 , vcov_wave)

coefs <- coeftest(m.0 , vcov_wave)[,1]
coefs <- coefs[-c(1, 13:16)]
ci.lo <- coeftest(m.0 , vcov_wave)[,1] - 1.96 * coeftest(m.0 , vcov_wave)[,2]
ci.hi <- coeftest(m.0 , vcov_wave)[,1] + 1.96 * coeftest(m.0 , vcov_wave)[,2]
ci.lo <- ci.lo[-c(1, 13:16)]
ci.hi <- ci.hi[-c(1, 13:16)]

points(coefs, 1:length(coefs)-.2, pch=19, cex=.8, col="black")
segments(ci.lo, 1:length(coefs)-.2, ci.hi, 1:length(coefs)-.2, col="black")

text(.20, 11, "old", col="red", cex=.8)
text(.20, 10.5, "new", col="black", cex=.8)

##############################################################################################





